In [1]:
%%file output/custom.css
body {
    font-size: 160%;
    font-family: Lato, Ariel, sans-serif !important;
}

div.slides {
    margin-top: -7%;
}

.left {
    width: 49%;
    float: left;
}

.right {
    width: 49%;
    float: right;
}

.centre {
    text-align: center;
}

h1 {
    text-align: center;
}

h2 {
    text-align: center;
}

div.prompt {
    display: none;
}

div.highlight {
    font-size: 85%; /* This Seems to give an approximate 80char line at 1024 */
}

div.output_html {
    font-size: 85%;
}

div.output_subarea {
    max-width: 100% !important;
}

div.output_png {
    text-align: center;
}

li {
    padding-top: 1em !important;
}

ul.logos {
    margin: 0px;
    padding: 0px;
    width: 100%;
}

ul.logos li {
    list-style: none;
    height:150px;

}

.output_stderr {
    display: none;
}

div.output_area .rendered_html img, div.output_area .rendered_html table {
    margin-left: auto;
    margin-right: auto;
}

.data-table {
    border-collapse: collapse;
}

.border-bottom {
    border-bottom: 1px solid #000;
}

.align-center {
    text-align: center;
}

.reveal i {
    font-size: inherit !important;
    font-style: italic !important;
}

.cite2c-dialog .twitter-typeahead {
    margin: 8px;
    width: 90%;
}

.cite2c-dialog .tt-input {
    width: 100%;
}

.cite2c-dialog .tt-dropdown-menu {
    background-color: #fff;
    border: 1px solid #ccc;
    border-radius: 4px;
}

.cite2c-dialog .tt-suggestion {
    padding: 4px;
}
.cite2c-dialog .tt-suggestion.tt-cursor {
    background-color: #eee;
}

.csl-entry {
    padding: 0.4em 0 !important;
}
Overwriting output/custom.css
In [2]:
import sys

sys.path.append('./python/')

from mayavi import mlab
from mayavi.tools.sources import vector_field, scalar_field
mlab.options.offscreen = True

import numpy as np

import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import astropy.units as u
import sunpy.map

import yt
import pysac.yt

%matplotlib inline

font_size = 20
pgf_with_latex = {
    "font.size": font_size,
    "axes.labelsize": font_size,  # LaTeX default is 10pt font.
    "legend.fontsize":font_size,
    "xtick.labelsize": font_size,
    "ytick.labelsize": font_size,
    "savefig.transparent": True
    }

matplotlib.rcParams.update(pgf_with_latex)
WARNING:traits.has_traits:DEPRECATED: traits.has_traits.wrapped_class, 'the 'implements' class advisor has been deprecated. Use the 'provides' class decorator.

Simulations of Magnetohydrodynamic Waves Driven by Photospheric Motions

Stuart J. Mumford

Supervisor: Robertus Erdélyi

Solar Physics & Space Plasma Research Centre (SP2RC), School of Mathematics and Statistics, The University of Sheffield

Problem - Coronal Heating


The solar atmosphere is too hot compared to understood heating mechanisms.

What are the unknown heating mechanisms?

  • Magnetic Reconnection
  • Magnetohydrodynamic (MHD) Waves
  • ...

Heating by MHD Waves

  • Waves generated in high-energy convection zone.
  • Waves propagate through the atmosphere, guided by magnetic field lines.
  • Waves then have to deposit energy in the atmosphere.

The mechanism by which the wave energy is converted into atmospheric heating in the atmosphere is unknown, but the properties of the wave behaviour is dependent on the wave mode.

How are these waves generated and what are their properties?

Numerical Simulations of the lower solar atmosphere

The Code

The code used is the Sheffield Advanced Code (SAC) (Shelyag, Fedun, and Erdélyi 2008).

SAC simulates perturbations on a static background, using a CD4 solver with hyper-diffusion and hyper-viscosity terms to stabalise the solution.

This makes it well-suited to solving wave perturbations on top of a highly stratified background such as the solar atmosphere.

The Model

To simulate wave excitation in the photosphere a numerical model of the solar atmosphere is needed.

Hydrostatic background from the VAL III model C (Vernazza, Avrett, and Loeser 1981):

In [3]:
import pysac.mhs_atmosphere as atm

#Read in the VAL3C model
empirical_data = atm.hs_atmosphere.read_VAL3c_MTW(MTW_file=False)


# Create a Z array at the interpolated resolution and interpolate.
ZZ = u.Quantity(np.linspace(empirical_data['Z'][0], empirical_data['Z'][-1], 128), unit=empirical_data['Z'].unit)
table = atm.hs_atmosphere.interpolate_atmosphere(empirical_data, ZZ, s=0)


# Create a figure and make space for the axes.
fig, ax = plt.subplots(gridspec_kw={'right':0.77, 'left':0.16, 'bottom':0.13}, figsize=(13,8))

# Shortcut all the Mm conversion.
Z = empirical_data['Z'].to(u.Mm)

lrho, = ax.plot(Z, empirical_data['rho'].quantity.si, 'x', color='blue')
lrho_i, = ax.plot(ZZ.to(u.Mm), table['rho'].quantity.si, color='blue')

ax2 = ax.twinx()
lp, = ax2.plot(Z, empirical_data['p'].to(u.Pa), 'x', color='green')
lp_i, = ax2.plot(ZZ.to(u.Mm), table['p'].to(u.Pa), color='green')


ax3 = ax.twinx()
ax3.spines["right"].set_position(("axes", 1.2))
lt, = ax3.plot(Z, empirical_data['T'].to(u.K), 'x', color='red')
lt_i, = ax3.plot(ZZ.to(u.Mm), table['T'].to(u.K), color='red')


# Set primary axes properties and labels
ax.semilogy()
ax.set_ylabel(r"Density [{}]".format(lrho._yorig.unit._repr_latex_()))
ax.set_xlabel(r"Height [{}]".format(lrho._xorig.unit._repr_latex_()))
ax.set_xlim(Z[0].value-0.1, Z[-1].value+0.1)


# Pressure Axis
ax2.semilogy()
ax2.set_ylabel(r"Pressure [{}]".format(lp._yorig.unit._repr_latex_()))


# Temp axis
ax3.set_ylabel(r"Temperature [{}]".format(lt._yorig.unit._repr_latex_()))

ax.set_xlim([-0.02,1.62])
ax3.set_ylim([3500,7500])
# Set the colours for the ticks and the labels.
ax.tick_params(axis='y', colors=lrho.get_color())
ax2.tick_params(axis='y', colors=lp.get_color())
ax3.tick_params(axis='y', colors=lt.get_color())

ax.yaxis.label.set_color(lrho.get_color())
ax2.yaxis.label.set_color(lp.get_color())
ax3.yaxis.label.set_color(lt.get_color())
fig
Out[3]:
In [4]:
from pysac.mhs_atmosphere.parameters.model_pars import mfe_setup as model_pars
import pysac.mhs_atmosphere as atm

# Cheeky Reset to Photosphere
model_pars['xyz'][4] = 0*u.Mm
#==============================================================================
# Build the MFE flux tube model using pysac
#==============================================================================
# model setup
scales, physical_constants = atm.units_const.get_parameters()
option_pars = atm.options.set_options(model_pars, False, l_gdf=True)
coords = atm.model_pars.get_coords(model_pars['Nxyz'], u.Quantity(model_pars['xyz']))

#interpolate the hs 1D profiles from empirical data source[s]
empirical_data = atm.hs_atmosphere.read_VAL3c_MTW(mu=physical_constants['mu'])
table = atm.hs_atmosphere.interpolate_atmosphere(empirical_data, coords['Zext'])
table['rho'] += table['rho'].min()*3.6

# calculate 1d hydrostatic balance from empirical density profile
# the hs pressure balance is enhanced by pressure equivalent to the
# residual mean coronal magnetic pressure contribution once the magnetic
# field has been applied
magp_meanz = np.ones(len(coords['Z'])) * u.one
magp_meanz *= model_pars['pBplus']**2/(2*physical_constants['mu0'])

# Make the vertical profile 3D
pressure_z, rho_z, Rgas_z = atm.hs_atmosphere.vertical_profile(coords['Z'], table, magp_meanz,
                                                               physical_constants, coords['dz'])

# Generate 3D coordinate arrays
x, y, z = u.Quantity(np.mgrid[coords['xmin']:coords['xmax']:1j*model_pars['Nxyz'][0],
                              coords['ymin']:coords['ymax']:1j*model_pars['Nxyz'][1],
                              coords['zmin']:coords['zmax']:1j*model_pars['Nxyz'][2]], unit=coords['xmin'].unit)

# Get default MFE flux tube parameters out of pysac
xi, yi, Si = atm.flux_tubes.get_flux_tubes(model_pars, coords, option_pars)

# Generate the 3D magnetic field
allmag = atm.flux_tubes.construct_magnetic_field(x, y, z, xi[0], yi[0], Si[0], model_pars, option_pars,
                                                 physical_constants, scales)
pressure_m, rho_m, Bx, By ,Bz, Btensx, Btensy = allmag

# local proc 3D mhs arrays
pressure, rho = atm.mhs_3D.mhs_3D_profile(z, pressure_z, rho_z, pressure_m, rho_m)
magp = (Bx**2 + By**2 + Bz**2) / (2.*physical_constants['mu0'])
energy = atm.mhs_3D.get_internal_energy(pressure, magp, physical_constants)

#### YT STUFF ####

magnetic_field_x = lambda field, data: data['mag_field_x']
yt.add_field(("gas","magnetic_field_x"), function=magnetic_field_x, units=yt.units.T.units)
magnetic_field_y = lambda field, data: data['mag_field_y']
yt.add_field(("gas","magnetic_field_y"), function=magnetic_field_y, units=yt.units.T.units)
magnetic_field_z = lambda field, data: data['mag_field_z']
yt.add_field(("gas","magnetic_field_z"), function=magnetic_field_z, units=yt.units.T.units)

# Add derived Fields
def magnetic_field_strength(field, data):
    return np.sqrt(data["mag_field_x"]**2 + data["mag_field_y"]**2 + data["mag_field_z"]**2)
yt.add_field(("gas","magnetic_field_strength"), function=magnetic_field_strength, units=yt.units.T.units)

#def alfven_speed(field, data):
#    return np.sqrt(2.*data['magnetic_pressure']/data['density'])
#yt.add_field(("gas","alfven_speed"), function=alfven_speed, units=(yt.units.m/yt.units.s).units)

bbox = u.Quantity([u.Quantity([coords['xmin'], coords['xmax']]),
                   u.Quantity([coords['ymin'], coords['ymax']]),
                   u.Quantity([coords['zmin'], coords['zmax']])]).to(u.m).value

# Now build a yt DataSet with the generated data:
data = {'mag_field_x':yt.YTQuantity.from_astropy(Bx.decompose()),
        'mag_field_y':yt.YTQuantity.from_astropy(By.decompose()),
        'mag_field_z':yt.YTQuantity.from_astropy(Bz.decompose()),
        'pressure': yt.YTQuantity.from_astropy(pressure.decompose()),
        'magnetic_pressure': yt.YTQuantity.from_astropy(magp.decompose()),
        'density': yt.YTQuantity.from_astropy(rho.decompose())}

ds = yt.load_uniform_grid(data, x.shape, length_unit='m', magnetic_unit='T',
                          mass_unit='kg', periodicity=[False]*3, bbox=bbox)
/home/stuart/GitHub/SWAT/pysac/pysac/mhs_atmosphere/mhs_model/flux_tubes.py:312: Warning: pbbal.max() = -0.111891563752 kg / (m s2)
  warnings.warn("pbbal.max() = {}".format(pbbal.max().decompose()), Warning)
yt : [INFO     ] 2015-12-02 13:56:45,420 Parameters: current_time              = 0.0
yt : [INFO     ] 2015-12-02 13:56:45,421 Parameters: domain_dimensions         = [129 129 128]
yt : [INFO     ] 2015-12-02 13:56:45,422 Parameters: domain_left_edge          = [-1000000. -1000000.        0.]
yt : [INFO     ] 2015-12-02 13:56:45,424 Parameters: domain_right_edge         = [ 1000000.   1000000.   1587786.3]
yt : [INFO     ] 2015-12-02 13:56:45,425 Parameters: cosmological_simulation   = 0.0

The Magnetic Flux Tube

The magnetic field model follows (Gent et al. 2013) and is constructed as a self-similar, cylindrically symmetric, expanding field.

In [5]:
slc = yt.SlicePlot(ds, 'x', 'density', origin='lower-center-domain', axes_unit='Mm')
slc.set_figure_size(10)
slc.set_cmap('all', 'viridis')
slc.set_font_size(20)
slc.set_zlim('all', 0, 2.5e-4)


seed_points = np.zeros([11,2]) + 1.52
seed_points[:,0] = np.linspace(-0.99, 0.95, seed_points.shape[0], endpoint=True)
slc.annotate_streamlines('mag_field_y', 'mag_field_z', field_color='magnetic_field_strength',
plot_args={'start_points':seed_points, 'density':15, 'cmap':'Blues', 'linewidth':2,
'norm': matplotlib.colors.LogNorm(*ds.all_data().quantities.extrema("magnetic_field_strength"))})
#force render
slc.save('/tmp/test.png')
#use the raw Figure for transparent bg
slc.plots['density'].figure
yt : [INFO     ] 2015-12-02 13:56:45,552 Loading field plugins.
yt : [INFO     ] 2015-12-02 13:56:45,553 Loaded angular_momentum (8 new fields)
yt : [INFO     ] 2015-12-02 13:56:45,554 Loaded astro (15 new fields)
yt : [INFO     ] 2015-12-02 13:56:45,556 Loaded cosmology (22 new fields)
yt : [INFO     ] 2015-12-02 13:56:45,558 Loaded fluid (63 new fields)
yt : [INFO     ] 2015-12-02 13:56:45,560 Loaded fluid_vector (95 new fields)
yt : [INFO     ] 2015-12-02 13:56:45,561 Loaded geometric (111 new fields)
yt : [INFO     ] 2015-12-02 13:56:45,562 Loaded local (115 new fields)
yt : [INFO     ] 2015-12-02 13:56:45,563 Loaded magnetic_field (121 new fields)
yt : [INFO     ] 2015-12-02 13:56:45,565 Loaded my_plugins (121 new fields)
yt : [INFO     ] 2015-12-02 13:56:45,566 Loaded species (123 new fields)
yt : [INFO     ] 2015-12-02 13:56:46,015 xlim = -1000000.000000 1000000.000000
yt : [INFO     ] 2015-12-02 13:56:46,016 ylim = 0.000000 1587786.300000
yt : [INFO     ] 2015-12-02 13:56:46,020 xlim = -1000000.000000 1000000.000000
yt : [INFO     ] 2015-12-02 13:56:46,020 ylim = 0.000000 1587786.300000
yt : [INFO     ] 2015-12-02 13:56:46,031 Making a fixed resolution buffer of (('gas', 'density')) 800 by 800
yt : [WARNING  ] 2015-12-02 13:56:46,101 Plot image for field ('gas', 'density') has both positive and negative values. Min = -0.000032, Max = 0.000267.
yt : [WARNING  ] 2015-12-02 13:56:46,102 Switching to symlog colorbar scaling unless linear scaling is specified later
yt : [INFO     ] 2015-12-02 13:56:48,090 Saving plot /tmp/test.png
Out[5]:

Exciting Waves in the Photosphere

In [6]:
x_range, y_range = [-300, 300]*u.arcsec, [-250, 250]*u.arcsec

plt.ioff()
fig = plt.figure(dpi=50, figsize=(11,8))
mm = sunpy.map.Map('/home/stuart/VivaData/gband_data/Gband_cospatial_cotemporal_00000.fits').submap(x_range, y_range)
mm = mm.submap([-440,440]*u.arcsec, [-440,440]*u.arcsec)
mm.plot_settings['cmap'] = 'gray'
mm.plot_settings['title'] = ''
im = mm.plot()
fig.savefig('./images/gband-plot.png', transparent=True)

The dynamic photosphere with embedded magnetic field provides many potential ways of driving MHD waves.

  • 'Buffeting' from convective motions
  • Global oscillations
  • Spiralling in downdrafts

Driving Waves in the Simulation Domain

$$ V(x,y,z) = F(x,y,z) \ e^{-\left(\frac{z^2}{\Delta z^2} + \frac{x^2}{\Delta x^2} + \frac{y^2}{\Delta y^2}\right)} \sin \left(2\pi \frac{t}{P}\right) $$

Identifying Waves from Broadband Drivers

Photospheric drivers excite multiple wave modes simulatenously.

How to quantify the relative strengths of the different modes from different drivers?


Assume uniform media:

  • Locally decompose perturbations into Fast, Slow and Alfvén modes.
  • Compare the percentage wave energy in each mode.
In [7]:
#Define tvtk notebook viewer
from IPython.core.display import Image
import subprocess
def mlab_view(scene, azimuth=153, elevation=62, distance=400, focalpoint=np.array([  25.,   63.,  60.]), aa=16):
    scene.anti_aliasing_frames = aa
    mlab.view(azimuth=azimuth, elevation=elevation, distance=distance, focalpoint=focalpoint)
    scene.save('offscreen.png', size=(750, 750))
    subprocess.call(["convert", "offscreen.png", "-transparent", "white", "offscreen.png"])
    return Image(filename='offscreen.png') 

Relationship between Mode and Velocity Perturbation

Assuming parallel propagation (along the magnetic field) $k_\parallel >> k_\perp$ it can be shown that the slow mode and fast mode perturb different components of velocity:

Fast mode: $$ \frac{|v_\parallel|}{|v_\perp|+|v_\parallel|} = \frac{\frac{k_\parallel}{k_\perp}|v_\perp|}{|v_\perp|+\frac{k_\parallel}{k_\perp}|v_\perp|}\\ = \frac{\frac{k_\parallel}{k_\perp}}{1+\frac{k_\parallel}{k_\perp}}\\ = \frac{1}{\frac{k_\perp}{k_\parallel}+1}\\ \ \\ |v_\parallel| >> |v_\perp| $$
Slow mode: $$ \frac{|v_\parallel|}{|v_\perp|+|v_\parallel|} = \frac{\frac{k_\perp}{k_\parallel}|v_\perp|}{|v_\perp|+\frac{k_\perp}{k_\parallel}|v_\perp|}\\ = \frac{\frac{k_\perp}{k_\parallel}}{1+\frac{k_\perp}{k_\parallel}}\\ = \frac{1}{\frac{k_\parallel}{k_\perp} + 1}\\ \ \\ |v_\parallel| << |v_\perp| $$

Decomposing Velocity Pertubation in SAC

In [8]:
ds = pysac.yt.SACGDFDataset('/home/stuart/VivaData/Slog_p240-0_A10_B005_00001.gdf')
yt : [WARNING  ] 2015-12-02 13:56:50,525 'field_units' was overridden by 'dataset_units/density_bg'
yt : [WARNING  ] 2015-12-02 13:56:50,526 'field_units' was overridden by 'dataset_units/density_pert'
yt : [WARNING  ] 2015-12-02 13:56:50,529 'field_units' was overridden by 'dataset_units/internal_energy_bg'
yt : [WARNING  ] 2015-12-02 13:56:50,531 'field_units' was overridden by 'dataset_units/internal_energy_pert'
yt : [WARNING  ] 2015-12-02 13:56:50,534 'field_units' was overridden by 'dataset_units/mag_field_x_bg'
yt : [WARNING  ] 2015-12-02 13:56:50,536 'field_units' was overridden by 'dataset_units/mag_field_x_pert'
yt : [WARNING  ] 2015-12-02 13:56:50,538 'field_units' was overridden by 'dataset_units/mag_field_y_bg'
yt : [WARNING  ] 2015-12-02 13:56:50,539 'field_units' was overridden by 'dataset_units/mag_field_y_pert'
yt : [WARNING  ] 2015-12-02 13:56:50,541 'field_units' was overridden by 'dataset_units/mag_field_z_bg'
yt : [WARNING  ] 2015-12-02 13:56:50,543 'field_units' was overridden by 'dataset_units/mag_field_z_pert'
yt : [WARNING  ] 2015-12-02 13:56:50,555 'field_units' was overridden by 'dataset_units/velocity_x'
yt : [WARNING  ] 2015-12-02 13:56:50,557 'field_units' was overridden by 'dataset_units/velocity_y'
yt : [WARNING  ] 2015-12-02 13:56:50,560 'field_units' was overridden by 'dataset_units/velocity_z'
yt : [INFO     ] 2015-12-02 13:56:50,597 Parameters: current_time              = 1.00339395144
yt : [INFO     ] 2015-12-02 13:56:50,599 Parameters: domain_dimensions         = [128 128 128]
yt : [INFO     ] 2015-12-02 13:56:50,600 Parameters: domain_left_edge          = [  781250.    781250.   3664122.1]
yt : [INFO     ] 2015-12-02 13:56:50,602 Parameters: domain_right_edge         = [  1.99218750e+08   1.99218750e+08   1.58778630e+08]
yt : [INFO     ] 2015-12-02 13:56:50,604 Parameters: cosmological_simulation   = 0.0
In [9]:
from tvtk.api import tvtk

#pysac imports
import pysac.yt
import pysac.analysis.tube3D.tvtk_tube_functions as ttf
import pysac.plot.mayavi_plotting_functions as mpf
from pysac.plot.mayavi_seed_streamlines import SeedStreamline
from pysac.plot.divergingcolourmaps import get_mayavi_colourmap
from pysac.analysis.tube3D.process_utils import get_yt_mlab

### Load in and Config ###

# loaded above
ds = pysac.yt.SACGDFDataset('/home/stuart/VivaData/Slog_p240-0_A10_B005_00001.gdf')
tube_r = 60

#if running this creates a persistant window just get it out of the way!
mlab.options.offscreen = True
fig = mlab.figure(bgcolor=(1, 1, 1))

cg = ds.index.grids[0]

#Slices
cube_slice = np.s_[:,:,:-5]
x_slice = np.s_[:,:,:,:-5]

#Define the size of the domain
xmax, ymax, zmax = np.array(cg['density'].to_ndarray()[cube_slice].shape) - 1
domain = {'xmax':xmax, 'ymax':ymax, 'zmax':zmax}

bfield, vfield = get_yt_mlab(ds, cube_slice, flux=False)

#Create a scalar field of the magntiude of the vector field
bmag = mlab.pipeline.extract_vector_norm(bfield, name="Field line Normals")
yt : [INFO     ] 2015-12-02 13:56:50,920 Loading field plugins.
yt : [INFO     ] 2015-12-02 13:56:50,922 Loaded angular_momentum (8 new fields)
yt : [INFO     ] 2015-12-02 13:56:50,924 Loaded astro (15 new fields)
yt : [INFO     ] 2015-12-02 13:56:50,925 Loaded cosmology (22 new fields)
yt : [INFO     ] 2015-12-02 13:56:50,927 Loaded fluid (63 new fields)
yt : [INFO     ] 2015-12-02 13:56:50,931 Loaded fluid_vector (95 new fields)
yt : [INFO     ] 2015-12-02 13:56:50,932 Loaded geometric (111 new fields)
yt : [INFO     ] 2015-12-02 13:56:50,934 Loaded local (115 new fields)
yt : [INFO     ] 2015-12-02 13:56:50,935 Loaded magnetic_field (120 new fields)
yt : [INFO     ] 2015-12-02 13:56:50,936 Loaded my_plugins (120 new fields)
yt : [INFO     ] 2015-12-02 13:56:50,937 Loaded species (122 new fields)
In [10]:
xc = domain['xmax']/2
yc = domain['ymax']/2
ti = 0
n = 100

surf_seeds = []
for theta in np.linspace(0, 2 * np.pi, n, endpoint=False):
    surf_seeds.append([tube_r * np.cos(theta + 0.5 * ti) + xc,
    tube_r * np.sin(theta + 0.5 * ti) + yc, domain['zmax']])

seeds = np.array(surf_seeds)
#Add axes:
axes, outline = mpf.add_axes(np.array(zip(ds.domain_left_edge,ds.domain_right_edge)).flatten()/1e8)

#Add seed points to plot:
seed_points = mlab.points3d(seeds[:,0], seeds[:,1], seeds[:,2],
color=(0.231, 0.298, 0.752), scale_mode='none',
scale_factor=1.5)

mlab_view(fig.scene)
/opt/miniconda/envs/thesis/lib/python2.7/site-packages/mayavi/tools/camera.py:288: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if focalpoint is not None and not focalpoint == 'auto':
Out[10]:
In [11]:
field_lines = SeedStreamline(seed_points = np.array(seeds))
bmag.add_child(field_lines)
field_lines.actor.mapper.scalar_visibility = False
field_lines.actor.property.color = (0,0,0)
field_lines.actor.property.line_width = 1.5

mlab_view(fig.scene)
Out[11]:
In [12]:
pd_seeds = ttf.make_circle_seeds(100, 60, **domain)
fieldlines, surface = ttf.create_flux_surface(bfield.outputs[0], pd_seeds)
surface.output.lines = None
flux_surface = mlab.pipeline.surface(surface.output)
flux_surface.actor.mapper.scalar_visibility = False
flux_surface.actor.property.color = (0.8,0.8,0.8)
#flux_surface.actor.property.line_width = 0

mlab_view(fig.scene)
Out[12]:
In [13]:
axes.visible = False
outline.visible = False
flux_surface.actor.property.edge_visibility = True
mlab_view(fig.scene, azimuth = 90, elevation = 75, distance=80, focalpoint=[63, 120, 110], aa=20)
Out[13]:
In [14]:
poly_norms = ttf.make_poly_norms(surface.output)
normvec = mlab.pipeline.glyph(poly_norms.output)
normvec.glyph.glyph_source.glyph_source = normvec.glyph.glyph_source.glyph_dict['arrow_source']
normvec.glyph.glyph.scale_mode = 'data_scaling_off'
normvec.glyph.glyph.color_mode = 'color_by_scale'
normvec.glyph.glyph.scale_factor = 5
normvec.glyph.glyph_source.glyph_position = 'tail'

mlab_view(fig.scene, azimuth=85, elevation=80, distance=50, focalpoint=[63, 120, 110], aa=20)
Out[14]:

MHD Waves excited by Different Photospheric Drivers

Chapter 4




Mumford, S. J., Fedun, V., Erdélyi, R. - The Astrophysical Journal - January 2015 - Volume 799, Issue 1
Generation of Magnetohydrodynamic Waves in Low Solar Atmospheric Flux Tubes by Photospheric Motions

In [15]:
from streamlines import Streamlines
#Use Equation 1 to calculate the vector field in a 2D plane to plot it.
time = np.linspace(0,60,480)
dt = time[1:] - time [:-1]
period = 240.

x = np.linspace(7812.5,1992187.5,128)
y = np.linspace(7812.5,1992187.5,128)

x_max = x.max()
y_max = y.max()

xc = 1.0e6
yc = 1.0e6

xn = x - xc
yn = y - yc

delta_x=0.1e6
delta_y=0.1e6

xx, yy = np.meshgrid(xn,yn)
exp_y = np.exp(-(yn**2.0/delta_y**2.0))
exp_x = np.exp(-(xn**2.0/delta_x**2.0))

exp_x2, exp_y2= np.meshgrid(exp_x,exp_y)
exp_xyz = exp_x2 * exp_y2


#==============================================================================
# Define Driver Equations and Parameters
#==============================================================================
#A is the amplitude, B is the spiral expansion factor
A = 10

#Tdamp defines the damping of the driver with time, Tdep is the ocillator
tdamp = lambda time1: 1.0 #*np.exp(-(time1/(period)))
tdep = lambda time1: np.sin((time1*2.0*np.pi)/period) * tdamp(time1)

#Define a peak index to use for scaling in the inital frame
max_ind = np.argmax(tdep(time) > 0.9998)

def log():
    B = 0.05
    phi = np.arctan2(1,B)
    theta = np.arctan2(yy,xx)

    uy = np.sin(theta + phi)
    ux =  np.cos(theta + phi)

    vx = lambda time1: (ux / np.sqrt(ux**2 + uy**2)) * exp_xyz * tdep(time1) * A
    vy = lambda time1: (uy / np.sqrt(ux**2 + uy**2)) * exp_xyz * tdep(time1) * A

    vv = np.sqrt(vx(time[max_ind])**2 + vy(time[max_ind])**2)

    return vx, vy, vv

def arch():
    B = 0.005
    r = np.sqrt(xx**2 + yy**2)

    vx = lambda time1: ( (B*1e6 * xx) / (xx**2 + yy**2) + yy/r ) * exp_xyz * tdep(time1) * A
    vy = lambda time1: ( (B*1e6 * yy) / (xx**2 + yy**2) - xx/r ) * exp_xyz * tdep(time1) * A

    vv = np.sqrt(vx(time[max_ind])**2 + vy(time[max_ind])**2)

    return vx, vy, vv

def uniform():
    #Uniform
    vx = lambda time1: A * (yy / np.sqrt(xx**2 + yy**2)) * exp_xyz * tdep(time1)
    vy = lambda time1: A * (-xx / np.sqrt(xx**2 + yy**2)) * exp_xyz * tdep(time1)
    vv = np.sqrt(vx(time[max_ind])**2 + vy(time[max_ind])**2)
    
    return vx, vy, vv

drivers = [log, arch, uniform]

#fig, axs = plt.subplots(1, 3, figsize=(18,9))

for driver_func in drivers:
    fig, ax = plt.subplots(figsize=(7,6), dpi=300)
    #============================================================================
    # Do the Plotting
    #============================================================================
    vx, vy, vv = driver_func()
    # Calculate Streamline
    slines = Streamlines(x,y,vx(time[max_ind]),vy(time[max_ind]),maxLen=7000,
                         x0=xc, y0=yc, direction='forwards')

    im = ax.imshow(vv, cmap='Blues', extent=[7812.5,x_max,7812.5,y_max])
    im.set_norm(matplotlib.colors.Normalize(vmin=0,vmax=A))
    #ax.hold()
    
    if driver_func != uniform:
        Sline, = ax.plot(slines.streamlines[0][0],slines.streamlines[0][1],color='red',linewidth=2, zorder=40)
    else:
        Sline = matplotlib.patches.Circle([1e6, 1e6], radius=.15e6, fill=False, color='red', linewidth=2, zorder=40)
        ax.add_artist(Sline)

    #Add colourbar
    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", size="5%", pad=0.2)
    cbar = plt.colorbar(im,cax)
    cbar.set_label(r"$|V|$ [ms$^{-1}$]")
    scalar = matplotlib.ticker.ScalarFormatter(useMathText=False,useOffset=False)
    scalar.set_powerlimits((-3,3))
    cbar.formatter = scalar
    cbar.ax.yaxis.get_offset_text().set_visible(True)
    cbar.update_ticks()
    #cbar.solids.set_rasterized(True)
    cbar.solids.set_edgecolor("face")

    #Add quiver plot overlay
    qu = ax.quiver(x, y, vx(time[max_ind]), vy(time[max_ind]), scale=25*A, color='k', zorder=20, linewidth=1)
    ax.axis([8.0e5,12.0e5,8.0e5,12.0e5])

    ax.xaxis.set_major_formatter(scalar)
    ax.yaxis.set_major_formatter(scalar)
    ax.xaxis.set_major_locator(matplotlib.ticker.MaxNLocator(5))
    ax.yaxis.set_major_locator(matplotlib.ticker.MaxNLocator(5))
    ax.xaxis.get_offset_text().set_visible(False)
    ax.yaxis.get_offset_text().set_visible(False)
    ax.set_xlabel("X [Mm]")
    ax.set_ylabel("Y [Mm]")

    fig.savefig('images/driver_{}.png'.format(driver_func.__name__), transparent=True)

Driving Waves in the Simulation Domain


$$ V_{(x,y,z)} = F_{(x,y,z)} \ e^{-\left(\frac{z^2}{\Delta z^2} + \frac{x^2}{\Delta x^2} + \frac{y^2}{\Delta y^2}\right)} \sin \left(2\pi \frac{t}{P}\right) $$

Uniform Driver

$$ F(x) = A \frac{y}{\sqrt{x^2 + y^2}},\\ F(y) = - A \frac{x}{\sqrt{x^2 + y^2}}, $$

Archmedian Spiral


$$ F(x) = A \frac{B_Ax}{x^2 + y^2} \frac{y}{\sqrt{x^2 + y^2}},\\ F(y) = - A \frac{B_Ay}{x^2 + y^2} \frac{x}{\sqrt{x^2 + y^2}}. $$

Logarithmic Spiral


$$ F(x) = A \frac{\cos(\theta + \phi)}{\sqrt{x^2 + y^2}},\\ F(y) = - A \frac{\sin(\theta + \phi)}{\sqrt{x^2 + y^2}},\\ $$ where, $\theta = \tan^{-1}\left(\frac{y}{x}\right),\ \phi = \tan^{-1}\left(\frac{1}{B_L}\right)$

Analysis and Results


Wave Flux

Calculate wave energy flux from (Leroy 1985). $$ \vec{F}_{\textbf{wave}} \equiv \widetilde{p}_k \vec{v} + \frac{1}{\mu_0} \left(\vec{B}_b \cdot \vec{\widetilde{B}}\right) \vec{v} - \frac{1}{\mu_0}\left(\vec{v} \cdot \vec{\widetilde{B}} \right) \vec{B}_b $$

Time-Distance Diagrams

In [16]:
import td_plotting_helpers as ph
import time_distance_plots as tdp
from matplotlib.image import NonUniformImage
import texfigure
ch4 = texfigure.Manager(None, number=4, base_path='/home/stuart/GitHub/Thesis/thesis/Chapter4/')

figsize = (17,8)

pvel_labels = {'par_label':r'$V_\parallel[$ ms$^{-1}$]', 
                'perp_label':r'$V_\perp$ [ms$^{-1}$]',
                'phi_label':r'$V_\phi$ [ms$^{-1}$]'}


pflux_labels = {'par_label':r'$F_\parallel / F^2$ ', 
                'perp_label':r'$F_\perp / F^2$',
                'phi_label':r'$F_\phi / F^2$'}

post_amp = "A10"
period = "p240"
tube_r = 'r30'
drivers = ['horiz', 'vert', 'Suni', 'Sarch', 'Slog']
exp_facs = [None, None, 'B0', 'B0005', 'B005']
captions = ['Horizontal', 'Vertical', 'Circular', 'Archemedian Spiral', 'Logarithmic Spiral']

figures = {}

for j, (driver, exp_fac, caption) in enumerate(zip(drivers, exp_facs, captions)):
    
    all_times, y, all_spoints = tdp.get_xy(ch4.data_dir, driver, period, post_amp, tube_r, exp_fac)
    data, beta_line = tdp.get_data(ch4.data_dir, driver, period, post_amp, tube_r, exp_fac)
    va_line, cs_line = tdp.get_speeds(ch4.data_dir, driver, period, post_amp, tube_r, exp_fac)

    Fdata, beta_line, avgs = tdp.get_flux(ch4.data_dir, driver, period, post_amp, tube_r, exp_fac)

    
    fd = lambda args: [a.T for a in args]
       
    ph.xxlim = -150
    
    fig, axes = plt.subplots(nrows=3, ncols=2, sharex=True, figsize=figsize)
    ph.triple_plot(axes[:,0], all_times, y, *fd(data), beta_line=1./beta_line,
                   levels=[1.,3.,5.,7.,10.,100.], manual_locations=False, **pvel_labels)
    
    
    
    ph.triple_plot(axes[:,1], all_times, y, *fd(Fdata), beta_line=1./beta_line,
                   levels=[1.,3.,5.,7.,10.,100.], manual_locations=False, cmap='PRGn', **pflux_labels)
    
    for ax in axes.flat:
        ph.add_phase_speeds(ax, all_times, y, va_line, cs_line, dx_scale=1e6, color='g')
        
    #Remove the top two x labels
    axes[0,0].set_xlabel('')
    axes[1,0].set_xlabel('')
    axes[0,1].set_xlabel('')
    axes[1,1].set_xlabel('')
    fig.tight_layout(h_pad=0.05)
    figures[driver] = fig
/opt/miniconda/envs/thesis/lib/python2.7/site-packages/matplotlib/__init__.py:1350: UserWarning:  This call to matplotlib.use() has no effect
because the backend has already been chosen;
matplotlib.use() must be called *before* pylab, matplotlib.pyplot,
or matplotlib.backends is imported for the first time.

  warnings.warn(_use_error_msg)

Vertical Driver

In [17]:
figures['vert']
Out[17]:

Horizontal Driver

In [18]:
figures['horiz']
Out[18]:

Logarithmic Spiral

In [19]:
figures['Slog']
Out[19]:

Comparison of Flux Averages

In [20]:
from flux_comparison import make_flux_bar_chart, get_averages

Favgs = get_averages(ch4.data_dir)

fig = make_flux_bar_chart((13,11), ch4.data_dir)
fig
Out[20]:
In [21]:
from sacconfig import SACConfig
cfg = SACConfig(cfg_file="./python/ch5_config.cfg")

BL = np.array([0.015, 0.05, 0.15, 0.45, 1.5])
In [22]:
fig, ax = plt.subplots(figsize=(9,2), dpi=600)
ax.plot(BL, np.ones(BL.size), 'x', markersize=10, mew=2)
ax.errorbar([0.15], [1], xerr=np.array([[-1*(0.15-1/(6.4-1.6)), 0.15+1/(6.4+1.6)]]).T, mew=2, elinewidth=2)
ax.semilogx()
ax.get_yaxis().set_visible(False)
ax.set_frame_on(False)
ax.get_xaxis().tick_bottom()
ax.xaxis.set_tick_params(width=2)
ax.xaxis.set_tick_params(width=2, which='minor')
ax.xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
ax.xaxis.set_ticks(BL)
xmin, xmax, ymin, ymax = ax.axis()
ax.add_artist(plt.Line2D((xmin, xmax), (ymin, ymin), color='black', linewidth=1.4))
l = ax.set_xlim([0.01, 2.0])
l = ax.set_xlabel(r'$B_L$', fontsize=18)

fig.tight_layout(h_pad=0.01)
fig.savefig('./images/ch5_space.png', transparent=True)

Effects of Logarithmic Spiral Expansion Factor

Chapter 5




Mumford, S. J. and Erdélyi, R. - Monthly Noticies of the Royal Astronomical Society - March 2015 - Volume 449 Issue 2.
Photospheric Logarithmic Velocity Spirals as MHD Wave Generation Mechanisms

Expansion Factor


$$ F(x) = A \frac{\cos(\theta + \phi)}{\sqrt{x^2 + y^2}},\\ F(y) = - A \frac{\sin(\theta + \phi)}{\sqrt{x^2 + y^2}},\\ $$ where, $\theta = \tan^{-1}\left(\frac{y}{x}\right),\ \phi = \tan^{-1}\left(\frac{1}{B_L}\right)$
Using the result of (Bonet et al. 2008):
$$ B_L^{-1} = 6.4 \pm 1.6 = B_L = 0.15 $$
In [23]:
beta = False
def add_triple_phase(ax, tube_r):
    ps = ph.get_phase_speeds(cfg, tube_r)
    for ax0 in ax:
        ph.add_phase_speeds(ax0, color='g', **ps)

bl_figures = {}
for j, bl in enumerate(BL):
    cfg.exp_fac = bl
    
    fig, ax = plt.subplots(nrows=3, ncols=2, sharex=True, figsize=figsize)
    
    kwargs = ph.get_single_velocity(cfg, 'r30', beta=beta)
    kwargs.update(pvel_labels)
    
    ph.triple_plot(ax[:,0], **kwargs)
    
    kwargs = ph.get_single_percentage_flux(cfg, 'r30', beta=beta)
    kwargs.update(pflux_labels)
    kwargs.update({'cmap': 'PRGn'})
    
    ph.triple_plot(ax[:,1], **kwargs)
    
    #Remove the top two x labels
    ax[0,0].set_xlabel('')
    ax[1,0].set_xlabel('')
    add_triple_phase(ax[:,0], 'r30')
    ax[0,1].set_xlabel('')
    ax[1,1].set_xlabel('')
    add_triple_phase(ax[:,1], 'r30')
    #add_all_bpert(ax, 'r30')
    fig.tight_layout(h_pad=0.05)
    bl_figures['B{}'.format(str(bl).replace('.',''))] = fig

More Time-Distance Diagrams

$B=0.015$

In [24]:
bl_figures['B0015']
Out[24]:

$B=0.15$

In [25]:
bl_figures['B015']
Out[25]:

$B=1.5$

In [26]:
bl_figures['B15']
Out[26]:

Average Wave Flux

In [27]:
import os
int_periods = np.floor(600./cfg.period)*180

def calc_avgs(tube_r):
    AvgsP = np.zeros([3,len(BL)])
    for i, bl in enumerate(BL):
        cfg.exp_fac = bl
        
        times = np.load(os.path.join(cfg.data_dir, 'Times_{}.npy'.format(cfg.get_identifier())))
        max_index = np.argmin(np.abs(int_periods - times))
        
        Fpar, Fperp, Fphi = map(np.load, ph.glob_files(cfg, tube_r, 'LineFlux*Fp*npy'))
        #Fpar, Fperp, Fphi = map(np.load, ph.glob_files(cfg, tube_r, '*vp*npy'))
        Fpar[np.abs(Fpar)<1e-5], Fperp[np.abs(Fperp)<1e-5], Fphi[np.abs(Fphi)<1e-5] = 0., 0., 0.
        Fpar, Fperp, Fphi = Fpar[:max_index,:], Fperp[:max_index,:], Fphi[:max_index,:]
        
        Ftot2 = (Fpar**2 + Fperp**2 + Fphi**2)
        Fpar2, Fperp2, Fphi2 = np.array([Fpar, Fperp, Fphi])**2
        FparP, FperpP, FphiP = (Fpar2/Ftot2)*100, (Fperp2/Ftot2)*100, (Fphi2/Ftot2)*100
        
        FparP = np.mean(np.ma.masked_array(FparP, np.isnan(FparP)))
        FperpP = np.mean(np.ma.masked_array(FperpP, np.isnan(FperpP)))
        FphiP = np.mean(np.ma.masked_array(FphiP, np.isnan(FphiP)))
        
        AvgsP[:, i] = FparP, FperpP, FphiP
    return AvgsP

figsize = (9.5,11)
fig, axs = plt.subplots(nrows=3, figsize=figsize, sharex=True)
titles = ["Flux Surface Radius $=156$ km", "Flux Surface Radius $=468$ km", "Flux Surface Radius $=936$ km"]
tubes = []
for i, ax in enumerate(axs):
    AvgsP = calc_avgs(cfg.tube_radii[i])
    tubes.append(AvgsP)
    ax.semilogx()
    ax.plot(BL, AvgsP[0], 'o', label=r"$F_\parallel^2$", mew=0, ms=10)
    ax.plot(BL, AvgsP[1], '_', label=r"$F_\perp^2$", mew=2, ms=10)
    ax.plot(BL, AvgsP[2], 'x', label=r"$F_\phi^2$", mew=2, ms=10)
    ax.set_ylabel("% Square Wave Flux")
    ax.set_title(titles[i])
    ax.xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
    ax.xaxis.set_ticks(BL)
    ax.set_xlim([0.01, 2.01])
    ax.set_ylim([0, 85])
    err = np.array([-1*(0.15-1/(6.4-1.6)), 0.15+1/(6.4+1.6)])
    ax.fill_betweenx(np.linspace(-5,105), err[0], err[1], alpha=0.3, color='green', linewidth=0)

axs[0].legend(loc=9)
axs[-1].set_xlabel("Logarithmic Spiral Expansion Factor ($B_L$)")

plt.tight_layout()
fig
Out[27]:

Effects of Period on MHD Wave Generation from a Logarithmic Spiral Driver

Chapter 6

Maintaining Consistent Energy Input


The amplitude is calculated using:
$$ A^2 \propto \frac{1}{P} $$
with $A=10$ $m/s$ and $P=240$ s as the reference point.

Period [seconds] Amplitude [ms$^{-1}$]
$30.0$ $20\sqrt{2}$
$60.0$ $20$
$90.0$ $20\sqrt{\frac{2}{3}}$
$120.0$ $10\sqrt{2}$
$150.0$ $4\sqrt{10}$
$180.0$ $\frac{20}{\sqrt{3}}$
$210.0$ $20\sqrt{\frac{2}{7}}$
$240.0$ $10$
$270.0$ $\frac{20}{3}\sqrt{2}$
$300.0$ $4\sqrt{5}$
In [28]:
from period_amps import periods, str_amps, sim_params
all_periods = sim_params[:10]
periods = periods[:10]

cfg = SACConfig(cfg_file='./python/ch6_config.cfg')
cfg.data_dir = '/home/stuart/GitHub/Thesis/thesis/Chapter6/Data/'
In [29]:
beta = False
cfg.exp_fac = 0.15
ph.xxlim = 600
tube_r = 'r30'

def add_all_bpert(ax, tube_r, N=4, levels=None):
    kwargs = ph.get_triple(cfg, beta=beta, single='bpert')
    x = kwargs['x_{}'.format(tube_r)]
    y = kwargs['y_{}'.format(tube_r)]
    par = kwargs['par_line_{}'.format(tube_r)].T[::-1, :]
    par[np.abs(par)<=1e-12] = 0
    perp = kwargs['perp_line_{}'.format(tube_r)].T[::-1, :]
    perp[np.abs(perp)<=1e-12] = 0
    phi = kwargs['phi_line_{}'.format(tube_r)].T[::-1, :]
    phi[np.abs(phi)<=1e-12] = 0
    ax[0].contour(x, y, par, N, colors='k', linewidths=np.linspace(0.5,1.5,N))
    ax[1].contour(x, y, perp, N, colors='k', linewidths=np.linspace(0.5,1.5,N))
    ax[2].contour(x, y, phi, N, colors='k', linewidths=np.linspace(0.5,1.5,N))	                   

def add_triple_phase(ax, tube_r):
    ps = ph.get_phase_speeds(cfg, tube_r)
    for ax0 in ax:
        ph.add_phase_speeds(ax0, color='g', **ps)

captions = {p: r"Period: ${}$ s amplitude: ".format(p) + a + r" ms$^{{-1}}$" for p, a in zip(periods, str_amps)[:10]}
#print(captions, file=sys.stderr)
width = 0.79

p_figures= {}
for i, paf in enumerate(all_periods):
    [setattr(cfg, f, getattr(paf, f)) for f in paf._fields]

    fig, ax = plt.subplots(nrows=3, ncols=2, figsize=(17,8), sharex=True)
    
    kwargs = ph.get_single_velocity(cfg, 'r30', beta=beta)
    kwargs.update(pvel_labels)
    
    ph.triple_plot(ax[:,0], **kwargs)
    
    kwargs = ph.get_single_percentage_flux(cfg, 'r30', beta=beta)
    kwargs.update(pflux_labels)
    kwargs.update({'cmap': 'PRGn'})
    
    ph.triple_plot(ax[:,1], **kwargs)
    
    #Remove the top two x labels
    ax[0,0].set_xlabel('')
    ax[1,0].set_xlabel('')
    add_triple_phase(ax[:,0], 'r30')
    ax[0,1].set_xlabel('')
    ax[1,1].set_xlabel('')
    add_triple_phase(ax[:,1], 'r30')
    #add_all_bpert(ax, 'r30')
    fig.tight_layout(h_pad=0.05)
    p_figures[str(paf.period)] = fig
/opt/miniconda/envs/thesis/lib/python2.7/site-packages/matplotlib/pyplot.py:516: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  max_open_warning, RuntimeWarning)

Even More Time-Distance Diagrams

$P=30$ $s$ $A=20\sqrt{2}$ $m/s$

In [30]:
p_figures['30.0']
Out[30]:

$P=150$ $s$ $A=4\sqrt{10}$ $m/s$

In [31]:
p_figures['150.0']
Out[31]:

$P=300$ $s$ $A=4\sqrt{5}$ $m/s$

In [32]:
p_figures['300.0']
Out[32]:

Period Comparison

In [33]:
from period_amps import periods, sim_params
sim_params = sim_params[:10]
periods = np.array(periods[:10])

fig, axs = plt.subplots(nrows=3, figsize=(10,11), sharex=True)
titles = ["Flux Surface Radius $=156$ km", "Flux Surface Radius $=468$ km", "Flux Surface Radius $=936$ km"]
tubes = []
for i, ax in enumerate(axs):
    AvgsP = ph.get_all_avgs(cfg, cfg.tube_radii[i], sim_params)
    tubes.append(AvgsP)
    ax.plot(periods, AvgsP[0], 'o', label=r"$F_\parallel^2$", mew=0, ms=12)
    ax.plot(periods, AvgsP[1], '_', label=r"$F_\perp^2$", mew=2, ms=12)
    ax.plot(periods, AvgsP[2], 'x', label=r"$F_\phi^2$", mew=2, ms=12)
    ax.set_ylabel("% Square Wave Flux")
    ax.set_title(titles[i])
    ax.xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
    ax.xaxis.set_ticks(periods)
    ax.set_xticklabels(["{}\n[{:n}]".format(p, 600//p) for p in periods])
    ax.set_ylim([10, 75])
    ax.set_xlim([25, 305])

axs[-1].set_xlabel("Period [s] \n [Number of periods averaged]")
axs[0].legend(loc=2)

#axs[0].legend(bbox_to_anchor=(1.06, 1.05))
plt.tight_layout(h_pad=0.1)
fig
Out[33]:

Conclusions

  • Different spiral driver profiles do not have a large impact on the excited waves.
  • Spirals with small expansion factors only excite $40-60$% Alfvén mode.
  • Spirals with large expansion factors exicte mostly fast mode.
  • Period only changes excited flux by $10-20$%

Future Work

Where we're going we don't need roads.

  • Vary the width of the driver FWHM.
  • Extend the background atmosphere using (Gent et al. 2013; Gent, Fedun, and Erdélyi 2014).
  • Decrease the magnetic field footpoint strength or change its size.
  • Multiple flux tubes for more realistic atmospheres.

Bibliography

Bonet, J. A., I. Márquez, J. Sánchez Almeida, I. Cabello, and V. Domingo. 2008. “Convectively Driven Vortex Flows in the Sun.” The Astrophysical Journal 687 (2): L131–34. doi:10.1086/593329.
Gent, F. A., V. Fedun, and R. Erdélyi. 2014. “Magnetohydrostatic Equilibrium. II. Three-Dimensional Multiple Open Magnetic Flux Tubes in the Stratified Solar Atmosphere.” The Astrophysical Journal 789 (1): 42. doi:10.1088/0004-637X/789/1/42.
Gent, F. A., V. Fedun, S. J. Mumford, and R. Erdelyi. 2013. “Magnetohydrostatic Equilibrium - I. Three-Dimensional Open Magnetic Flux Tube in the Stratified Solar Atmosphere.” Monthly Notices of the Royal Astronomical Society 435 (1): 689–97. doi:10.1093/mnras/stt1328.
Leroy, Bernard. 1985. “On the Derivation of the Energy Flux of Linear Magnetohydrodynamic Waves.” Geophysical & Astrophysical Fluid Dynamics 32 (2): 123–33. doi:10.1080/03091928508208781.
Shelyag, S., V. Fedun, and R. Erdélyi. 2008. “Magnetohydrodynamic Code for Gravitationally-Stratified Media.” Astronomy and Astrophysics 486 (2): 655–62. doi:10.1051/0004-6361:200809800.
Vernazza, J. E., E. H. Avrett, and R. Loeser. 1981. “Structure of the Solar Chromosphere. III - Models of the EUV Brightness Components of the Quiet-Sun.” The Astrophysical Journal Supplement Series 45 (April): 635. doi:10.1086/190731.

My Publications


Mumford, S. J. and Erdélyi, R. - Monthly Noticies of the Royal Astronomical Society - March 2015 - Volume 449 Issue 2.
Photospheric Logarithmic Velocity Spirals as MHD Wave Generation Mechanisms

Mumford, S. J., Fedun, V., Erdélyi, R. - The Astrophysical Journal - January 2015 - Volume 799, Issue 1
Generation of Magnetohydrodynamic Waves in Low Solar Atmospheric Flux Tubes by Photospheric Motions

The SunPy Community, Mumford, S. J., Christe, S., Pérez-Suárez, D., et. al - Computational Science and Discovery - January 2015 - Volume 8 Issue 1.
SunPy: Python for Solar Physics

Freij N., Scullion E. M., Nelson C. J., Mumford S. J., Wedemeyer S., and Erdélyi R. - The Astrophysical Journal - July 2014 - Volume 791, Issue 1, p.61
The Detection of Upwardly Propagating Waves Channeling Energy from the Chromosphere to the Low Corona

Gent, F. A., Fedun, V., Mumford, S. J., Erdélyi, R. - Monthly Notices of the Royal Astronomical Society - October 2013 - Volume 435, Issue 1, p.689-697
Magnetohydrostatic equilibrium - I. Three-dimensional open magnetic flux tube in the stratified solar atmosphere

Nelson, C. J., Doyle, J. G., Erdélyi, R., Huang, Z., Madjarska, M. S., Mathioudakis, M., Mumford, S. J., Reardon, K - Solar Physics - April 2013 - Volume 283, Issue 2, p.307-323.
Statistical Analysis of Small Ellerman Bomb Events